#!/Pub/Apps/bin/Rscript
suppressPackageStartupMessages(library(argparse))

wrp <- "转化pySCENIC分析之后的结果loom文件，提取其中的regulons以及regulonsAUC，转化为rds，保存在loom结果文件的同级目录下，为后续在R中可视化做准备，需要结果loom文件路径。"
parser <- ArgumentParser(
    prog='seurat2loom_exprMat.r',
    description = wrp,
    formatter_class = "argparse.RawTextHelpFormatter"
    # usage = "seurat2loom_exprMat.r -h; seurat2loom_exprMat.r seurat_obj od"
)

parser$add_argument("SCENIC_loom_path", nargs = 1, help = "SCEINC分析结果路径。[.loom file]")

args_in_process <- parser$parse_args()
cli::cli_inform(">> 1.1 Load pySCEINC loom obj [{Sys.time()}]\n")

od <- dirname(args_in_process$SCENIC_loom_path)
loom_path <- args_in_process$SCENIC_loom_path

suppressPackageStartupMessages(library(Seurat))
suppressPackageStartupMessages(library(SCopeLoomR))
suppressPackageStartupMessages(library(SCENIC))
suppressPackageStartupMessages(library(AUCell))
suppressPackageStartupMessages(library(tidyverse))

Sys.setenv(HDF5_USE_FILE_LOCKING = "FALSE")

main <- \(){
    loom <- SCopeLoomR::open_loom(loom_path)
    regulons_incidMat <- get_regulons(loom, column.attr.name = "Regulons")
    regulons <- regulonsToGeneLists(regulons_incidMat)
    regulonAUC <- get_regulons_AUC(loom, column.attr.name = "RegulonsAUC")
    regulonAucThresholds <- get_regulon_thresholds(loom)
    # tail(regulonAucThresholds[order(as.numeric(names(regulonAucThresholds)))])

    embeddings <- get_embeddings(loom)
    close_loom(loom)

    AUCellThresholds_Info <- AUCell_exploreThresholds(
        cellsAUC = regulonAUC,
        thrP = 0.01,
        nCores = 10,
        smallestPopPercent = 0.25,
        plotHist = TRUE,
        densAdjust = 2,
        assignCells = FALSE,
        nBreaks = 100,
        verbose = TRUE
    )
    # class(AUCellThresholds_Info)

    map(names(AUCellThresholds_Info), ~ {
        tibble(regulon = .x, threshold = AUCellThresholds_Info[[.x]]$aucThr$selected)
    }) %>% bind_rows() -> AUCellThresholds_Info


    saveRDS(list("regulons" = regulons, "regulonAUC" = regulonAUC, "AUCellThresholds_Info" = AUCellThresholds_Info),
        file = paste0(od, "/", "pyscenic_output.rds")
    )
    # 一些可视化参考 http://www.360doc.com/content/22/1104/09/74726315_1054515292.shtml
}

try(main())
